﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace TopRelax2
{
    public class SORelax
    {
        //код для простой итерации
        public double[,] u; //искомое решение
        public double[,] b;
        public int n; //число разбиений по икс
        public int m; //число разбиений по игрек
        public double h; // шаг по икс
        public double k; // шаг по игрек
        
        public double w; // параметр метода

        private double di;

        public SORelax()
        {
            n = m = 0;
            h = k = 0;
            setDi();
            u = null;
            b = null;
            
            w = 0;
        }
        public SORelax(int n, double h, int m, double k, double[,] b, double[,] u)
        {
            this.n = n;
            this.m = m;
            this.h = h;
            this.k = k;
            setDi();
            this.u = new double[2 * n + 1, 2 * m + 1];
            this.b = new double[2 * n + 1, 2 * m + 1];
            for (int i = 0; i <= 2 * n; i++)
            {
                for (int j = 0; j <= 2 * m; j++)
                {
                    this.u[i, j] = u[i, j];
                    this.b[i, j] = b[i, j];
                }
            }
            
            setW();

        }
        private void setDi()
        {
            di = -2.0 * (1.0 / (h * h) + 1.0 / (k * k));
        }
        private void setW()
        {
            int rank = (2 * n - 1) * (m - 1) + (n - 1) * m;
            double[,] matr = new double[rank, rank];
            double d1 = 1.0 / (h * h);
            double dn = 1.0 / (k * k);

            //заполнение матрицы
            int lineCounter = 0;
            for (int j = 1; j < m; j++)
            {
                for (int i = 1; i < 2 * n; i++)
                {
                    matr[lineCounter, lineCounter] = 0;
                    if (i != 1) matr[lineCounter, lineCounter - 1] = d1/di;
                    if (i != 2 * n - 1) matr[lineCounter, lineCounter + 1] = d1/di;
                    if (j != 1) matr[lineCounter, lineCounter - 2 * n + 1] = dn/di;
                    if (j != m - 1)
                    {
                        matr[lineCounter, lineCounter + 2 * n - 1] = dn/di;
                    }
                    else
                    {

                        if (i < n - 1)
                        {
                            matr[lineCounter, lineCounter + 2 * n - 1] = dn/di;
                        }
                    }
                    lineCounter++;
                }
            }
            for (int j = m; j < 2 * m; j++)
            {
                for (int i = 1; i < n; i++)
                {
                    matr[lineCounter, lineCounter] = 0;
                    if (i != 1) matr[lineCounter, lineCounter - 1] = d1/di;
                    if (i != n - 1) matr[lineCounter, lineCounter + 1] = d1/di;
                    if (j != 2 * m - 1) matr[lineCounter, lineCounter + n - 1] = dn/di;
                    if (j != m)
                    {
                        matr[lineCounter, lineCounter - n + 1] = dn/di;
                    }
                    else
                    {
                        matr[lineCounter, lineCounter - 2 * n + 1] = dn/di;
                    }
                    lineCounter++;
                }
            }

            //поиск собственных чисел
            double[] wr =new double[rank];
            double[] wi;
            double[,] wl,wu;
            alglib.rmatrixevd(matr, rank, 0, out wr, out wi,out wl,out wu);

            //поиск максимума
            int Imax = 0;
            for (int i = 0; i < rank; i++)
            {
                if (Math.Abs(wr[Imax]) < Math.Abs(wr[i])) Imax = i;   
            }

            //формула оптимального параметра
            w = 2.0 / (1 + Math.Sqrt(1 - wr[Imax] * wr[Imax]));
        }

        
        private double findMaxAbs(double[,] d1, double[,] d2, out int Imax, out int Jmax)
        {
            Imax = 1;
            Jmax = 1;
            for (int i = 1; i < 2 * n - 1; i++)
                for (int j = 1; j < 2 * m - 1; j++)
                    if (Math.Abs(d1[i, j] - d2[i, j]) > Math.Abs(d1[Imax, Jmax] - d2[Imax, Jmax]))
                    {
                        Imax = i;
                        Jmax = j;
                    }

            return Math.Abs(d1[Imax, Jmax] - d2[Imax, Jmax]);
        }

        public int method(double e, out double ereal, int count, int Imax, int Jmax)
        {

            int s = 0;
            double[,] oldu = new double[2 * n + 1, 2 * m + 1];

            double dl = 10;

            double d1 = 1.0 / (h * h);
            double dn = 1.0 / (k * k);

            while (dl > e && s < count)
            {
                s++;

                
                    for (int p = 0; p <= 2 * n; p++)
                        for (int j = 0; j <= 2 * m; j++)
                            oldu[p, j] = u[p, j];

                    //начало нового лгоритма
                    for (int j = 1; j < m; j++)
                    {
                        for (int p = 1; p < 2 * n; p++)
                        {
                            u[p, j] = (w * (-d1 * (u[p - 1, j] + oldu[p + 1, j]) - dn * (u[p, j - 1] + oldu[p, j + 1]) + b[p, j])) / di + (1 - w) * oldu[p, j];
                        }
                    }
                    for (int j = m; j < 2 * m; j++)
                    {
                        for (int p = 1; p < n; p++)
                        {
                            u[p, j] = (w * (-d1 * (u[p - 1, j] + oldu[p + 1, j]) - dn * (u[p, j - 1] + oldu[p, j + 1]) + b[p, j])) / di + (1 - w) * oldu[p, j];
                        }
                    }
                
                dl = findMaxAbs(u, oldu, out Imax, out Jmax);
            }
            ereal = dl;
            return s;
        }


    }
}
